#include "Lagrange.h"

#include <cassert>
#include <fstream>
#include <iostream>
#include <sstream>

#include "../../utils/MatrixUtils.h"

namespace NumericalAnalysis
{
    // Calculates the Lagrange polynomial li(x)
    double CalculateLagrangePolynomial(const Utils::Matriz& pairs, const double x, const uint32_t i)
    {
    	assert(pairs.mRows > 0 && pairs.mColumns == 2 
			&& "CalculateLagrangePolynomial: This is not a matrix of pairs.");

    	// Create the result.
    	double result = 1.0;
    	
    	// Update the result.
		// Multiply[j = 0; j != i; n] (x - Xj) / (Xi - Xj)
    	for(uint32_t j = 0; j < pairs.mRows; ++j)
    	{
    		if(j != i)
    		{
    			const double value1 = pairs.mData[j * pairs.mColumns]; // Xj
    			const double value2 = pairs.mData[i * pairs.mColumns]; // Xi
    			result *= (x - value1) / (value2 - value1); 
    		}
    	}

    	return result;
    }

    // Returns the pairs of (xi, Ln(xi)), where Ln 
    // is the Lagrange polynomial of interpolation.
    Utils::Matriz LagrangeInterpolation(const std::string& pairsFilePath, const uint32_t n, 
                                        const std::string& valuesFilePath, const uint32_t m)
    {
    	// Get the pairs: [xi, yi]
    	Utils::Matriz pairs = Utils::GetMatrixFromFile(pairsFilePath, n, 2);

    	// Get the values.
    	Utils::Matriz values = Utils::GetMatrixFromFile(valuesFilePath, m, 1);

    	// Create the pairs to return.
    	Utils::Matriz returnPairs(m, 2);
    	double *returnPairsData = returnPairs.mData;
    	// Calculate each pair.
    	for(uint32_t i = 0; i < values.mRows; ++i)
    	{
    		// Compute li(xi)
    		double newY = 0.0;
    		const double iValued = values.mData[i]; 
    	
    		for(uint32_t j = 0; j < pairs.mRows; ++j)
    		{
				// Ln(xi) = yj * lj(xj)
    			newY += pairs.mData[j * pairs.mColumns + 1] * CalculateLagrangePolynomial(pairs, iValued, j);
    		}

    		// Add the new pair.
    		returnPairsData[0] = iValued; // xi
    		returnPairsData[1] = newY; // Ln(xi)
    		returnPairsData += returnPairs.mColumns;
    	}

    	return returnPairs;
    }
}
